library(dplyr)
library(factoextra)
library(FactoMineR)
library(geosphere)
library(gridExtra)
library(leaflet)
library(lubridate)
library(mice)
library(patchwork)
library(purrr)
library(RColorBrewer)
library(readr)
library(sf)
library(tibble)
library(tidyr)

1 Mapa para visualizar conclusiones

trafico_final <- read_csv('/Users/pablogandia/Desktop/trafico_cluster.csv')
contaminantes <- read_csv("/Users/pablogandia/Desktop/Analisis_aire/data/raw/geoespacial/estaciones_valencia.csv") %>%
  select(Estación, Latitud, Longitud, Altitud)

zonas <- c("València - Av. França",
           "València - Molí del Sol", "València - Pista de Silla",
           "València - Politècnic", "València Port llit antic Túria",
           "València Port Moll Trans. Ponent")

contaminantes <- contaminantes[contaminantes$Estación %in% zonas, ]

datos_vias <- read_csv('/Users/pablogandia/Desktop/trafico_horario_clean.csv') %>%
  select(IdPM, Hora, Dia,Mes,IntensidadMediana,VelocidadMedia,OcupacionMedia,OcupacionMaxima,DevIntensidad)

pal <- colorFactor(palette = brewer.pal(3, "Set1"),
                   domain = trafico_final$tipo_via)

m <- leaflet(trafico_final) %>%
  addTiles() %>%
  setView(lng = mean(trafico_final$Lon), lat = mean(trafico_final$Lat), zoom = 11) %>%
  addCircleMarkers(
    lng = ~Lon, lat = ~Lat,
    radius = 6,
    color = ~pal(tipo_via),
    fillOpacity = 0.8,
    stroke = FALSE,
    popup = ~paste("Zona:", idpm, "<br>Tipo de vía:", tipo_via)
  ) %>%
  addLegend(
    position = "bottomright",
    pal = pal,
    values = ~tipo_via,
    title = "Tipos de vía",
    opacity = 1
  ) %>%
  addCircles(
    data = contaminantes,
    lng = ~Longitud,
    lat = ~Latitud,
    radius = 1000,
    color = "#ff9999",
    fillOpacity = 0.2,
    stroke = TRUE
  ) %>%
  addCircleMarkers(
    data = contaminantes,
    lng = ~Longitud,
    lat = ~Latitud,
    radius = 8,
    color = "black",
    fillOpacity = 1,
    popup = ~Estación
  )

m

2 Avenida de Francia

datos <- read_csv("/Users/pablogandia/Desktop/Trafico final 1000m/València_-_Av._França.csv")

datos <- datos %>% distinct(Mes, Dia, Hora, .keep_all = TRUE)

datos
vars_pca <- c(
  "NO", "NO2", "NOx", "O3", "PM10", "PM2.5", "SO2",
  "viento_Norte", "viento_Sur", "viento_Este", "viento_Oeste",
  "viento_Noreste", "viento_Sureste", "viento_Suroeste", "viento_Noroeste"
)

datos_pca <- datos %>%
  select(any_of(vars_pca)) %>%
  drop_na() %>%
  scale()

res_pca <- PCA(datos_pca, scale.unit = TRUE, graph = FALSE, ncp = 10)

eig.val <- get_eigenvalue(res_pca)
VPmedio <- 100 / nrow(eig.val)

fviz_eig(res_pca, addlabels = TRUE) +
  geom_hline(yintercept = VPmedio, linetype = 2, color = "red")

fviz_pca_var(res_pca,
             axes = c(1, 2),           
             col.var = "contrib",      
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE,             
             title = "Loadings: Componentes 1 y 2")

fviz_contrib(res_pca, choice = "var", axes = 1, top = Inf, 
             title = "Contribución de variables a la Componente 1")

fviz_contrib(res_pca, choice = "var", axes = 2, top = Inf, 
             title = "Contribución de variables a la Componente 2")

Vemos una relacion en la primera componente del viento suroeste con los óxidos de nitrogeno, también el viento del sur y del noroeste contribuyen a los PM en la primera y segunda componente, es decir, están relacionados con el valor de estos contaminantes.

Si visualizamos el mapa nos percatamos que justo estos son (principalmente noroeste y suroeste) los sentidos donde circulan vías de alta intensidad, por lo que el viento podría estar arrastrando las emisiones de los vehiculos, empeorando la calidad del aire.

vars_pca2 <- c(
  "NO", "NO2", "NOx", "O3", "PM10", "PM2.5", "SO2",
  "IntensidadGlobal", "VelocidadGlobal", "OcupacionGlobal")

datos_pca <- datos %>%
  select(any_of(vars_pca2)) %>%
  drop_na() %>%
  scale()

res_pca <- PCA(datos_pca, scale.unit = TRUE, graph = FALSE, ncp = 10)

eig.val <- get_eigenvalue(res_pca)
VPmedio <- 100 / nrow(eig.val)

fviz_eig(res_pca, addlabels = TRUE) +
  geom_hline(yintercept = VPmedio, linetype = 2, color = "red")

fviz_pca_var(res_pca,
             axes = c(1, 2),           
             col.var = "contrib",      
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE,             
             title = "Loadings: Componentes 1 y 2")

fviz_contrib(res_pca, choice = "var", axes = 1, top = Inf, 
             title = "Contribución de variables a la Componente 1")

fviz_contrib(res_pca, choice = "var", axes = 2, top = Inf, 
             title = "Contribución de variables a la Componente 2")

Contrariamente a lo que planteamos en base a la dirección del viento, el PCA indica que no existe relación lineal entre tráfico y contaminantes, y esque como habiamos adelantado en analisis previos, el trafico está tan relacionado entre si en el computo general que se junta en una sola componente, y no se relaciona con contaminación debido al granulo y contexto utilizado (los contaminantes en su mayoría forman una sola componente ya que como se ha visto están intensamente relacionados).

vars_pca2 <- c(
  "IntensidadGlobal", "VelocidadGlobal", "OcupacionGlobal",
  "Intensidad_alto", "Intensidad_bajo", "Intensidad_congestionado",
  "Velocidad_alto", "Velocidad_bajo", "Velocidad_congestionado",
  "Ocupacion_alto", "Ocupacion_bajo", "Ocupacion_congestionado")

datos_pca <- datos %>%
  select(any_of(vars_pca2)) %>%
  drop_na() %>%
  scale()

res_pca <- PCA(datos_pca, scale.unit = TRUE, graph = FALSE, ncp = 10)

eig.val <- get_eigenvalue(res_pca)
VPmedio <- 100 / nrow(eig.val)

fviz_eig(res_pca, addlabels = TRUE) +
  geom_hline(yintercept = VPmedio, linetype = 2, color = "red")

fviz_pca_var(res_pca,
             axes = c(1, 2),           
             col.var = "contrib",      
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE,             
             title = "Loadings: Componentes 1 y 2")

fviz_contrib(res_pca, choice = "var", axes = 1, top = Inf, 
             title = "Contribución de variables a la Componente 1")

fviz_contrib(res_pca, choice = "var", axes = 2, top = Inf, 
             title = "Contribución de variables a la Componente 2")

Es inutil utilizar otras variables de trafico generadas por nosotros, ya que como habiamos supuesto antes, están intensamente realacionadas todas, de hecho la primera componente explica el 92%, lo cual es ridiculamente alto.

3 Politécnico

datos <- read_csv("/Users/pablogandia/Desktop/Trafico final 1000m/València_-_Politècnic.csv",show_col_types = FALSE)

datos <- datos %>% distinct(Mes, Dia, Hora, .keep_all = TRUE)
vars_pca <- c(
  "NO", "NO2", "NOx", "O3", "PM10", "PM2.5", "SO2",
  "viento_Norte", "viento_Sur", "viento_Este", "viento_Oeste",
  "viento_Noreste", "viento_Sureste", "viento_Suroeste", "viento_Noroeste"
)

datos_pca <- datos %>%
  select(any_of(vars_pca)) %>%
  drop_na() %>%
  scale()

res_pca <- PCA(datos_pca, scale.unit = TRUE, graph = FALSE, ncp = 10)

eig.val <- get_eigenvalue(res_pca)
VPmedio <- 100 / nrow(eig.val)

fviz_eig(res_pca, addlabels = TRUE) +
  geom_hline(yintercept = VPmedio, linetype = 2, color = "red")

fviz_pca_var(res_pca,
             axes = c(1, 2),           
             col.var = "contrib",      
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE,             
             title = "Loadings: Componentes 1 y 2")

fviz_contrib(res_pca, choice = "var", axes = 1, top = Inf, 
             title = "Contribución de variables a la Componente 1")

fviz_contrib(res_pca, choice = "var", axes = 2, top = Inf, 
             title = "Contribución de variables a la Componente 2")

En este caso la dirección del viento que más se relaciona con los principales contaminantes, los óxidos de nitrógeno y los PM es el suroeste, que si nos fijamos en el mapa es justo la dirección que arrastraria las emisiones de la principal avenida cercana (avenida de los naranjos) hacia la estacion, que de nuevo implicaría una fuerte relación entre los contaminantes y el tráfico.

vars_pca2 <- c(
  "NO", "NO2", "NOx", "O3", "PM10", "PM2.5", "SO2",
  "IntensidadGlobal", "VelocidadGlobal", "OcupacionGlobal")

datos_pca <- datos %>%
  select(any_of(vars_pca2)) %>%
  drop_na() %>%
  scale()

res_pca <- PCA(datos_pca, scale.unit = TRUE, graph = FALSE, ncp = 10)

eig.val <- get_eigenvalue(res_pca)
VPmedio <- 100 / nrow(eig.val)

fviz_eig(res_pca, addlabels = TRUE) +
  geom_hline(yintercept = VPmedio, linetype = 2, color = "red")

fviz_pca_var(res_pca,
             axes = c(1, 2),           
             col.var = "contrib",      
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE,             
             title = "Loadings: Componentes 1 y 2")

fviz_contrib(res_pca, choice = "var", axes = 1, top = Inf, 
             title = "Contribución de variables a la Componente 1")

fviz_contrib(res_pca, choice = "var", axes = 2, top = Inf, 
             title = "Contribución de variables a la Componente 2")

No obstante los datos vuelven a mostrar una relación nula entre tráfico y contaminación, al menos bajo este paradigma.

4 Pista de Silla

datos <- read_csv("/Users/pablogandia/Desktop/Trafico final 1000m/València_-_Pista_de_Silla.csv",show_col_types = FALSE)

datos <- datos %>% distinct(Mes, Dia, Hora, .keep_all = TRUE)
vars_pca2 <- c(
  "NO", "NO2", "NOx", "O3", "PM10", "PM2.5", "SO2",
  "IntensidadGlobal", "VelocidadGlobal", "OcupacionGlobal")

datos_pca <- datos %>%
  select(any_of(vars_pca2)) %>%
  drop_na() %>%
  scale()

res_pca <- PCA(datos_pca, scale.unit = TRUE, graph = FALSE, ncp = 10)

eig.val <- get_eigenvalue(res_pca)
VPmedio <- 100 / nrow(eig.val)

fviz_eig(res_pca, addlabels = TRUE) +
  geom_hline(yintercept = VPmedio, linetype = 2, color = "red")

fviz_pca_var(res_pca,
             axes = c(1, 2),           
             col.var = "contrib",      
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE,             
             title = "Loadings: Componentes 1 y 2")

fviz_contrib(res_pca, choice = "var", axes = 1, top = Inf, 
             title = "Contribución de variables a la Componente 1")

fviz_contrib(res_pca, choice = "var", axes = 2, top = Inf, 
             title = "Contribución de variables a la Componente 2")

En la pista de silla, la principal salida y entrada de valencia, que presenta grandes volumenes de tráfico todos los dias ¡¡TAMPOCO!! muestra relación entre tráfico y contaminación, lo que nos hace concluir que el problema de este estudio a sido probablemente el enfoque. Quizas para explicar el trafico en una zona no basta con tomar medias horarias, ya que como hemos visto a lo largo del objetivo, las tendencias de tráfico son muy extensas y cambiantes, por lo que jamás encontraremos relaciones estáticas.

5 Conclusiones

Aunque algunos factores como la localización de la estación de contaminación respecto a las vías de trafico más cercanas, la dirección del viento y el propio sentido común nos indica que existe una fuerte relación entre el tráfico y la contaminación, el contexto bajo el que se pone el estudio hace que no se obtengan relaciones concluyentes, lo que refuerza lo anteriormente visto de como el tráfico varía mucho según la zona, la hora del dia, el dia de la semana…